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1 INTRODUCTION 



The Ro ssby w ave instability (RWI, iLovelace at al. 1 1 19991 : 
iLi at al.l I2OO0I ) is tha thin-disc analog of the Papaloizou- 
Pring l a inst ability in thick tori (jPapaloizou fc Pringle|[l983 . 
1 19851 . 119871 ). Tha RWI is a shear instability associated 
with extrema in the potential vorticity profile of a disc. 
Its occurrence in radially str uctured protq planetary discs, 
leading to vortex formation (|Li et al.ll2001^ ■ ca n transport 
angu l ar m omentum, assist planet formati on (iLvra et al.l 
2008, 2009) and affect planetary migration (|Li et aLlbOOa 
Lin fc Papaloizou 2OIQ1 ). 



The linear Rossby wave instability (RWI) in global, three-dimensional (3D) poly- 
tropic discs is revisited with a much simpler numerical method than that previously 
employed by the author. The governing partial differential equation is solved with 
finite-differences in the radial direction and spectral collocation in the vertical direc- 
tion. RWI modes are calculated subject to different upper disc boundary conditions. 
These include free surface, solid boundaries and variable vertical domain size. Bound- 
ary conditions that oppose vertical motion increases the instability growth rate by 
a few per cent. The magnitude of vertical flow throughout the fluid column can be 
affected but the overall flow pattern is qualitatively unchanged. Numerical results 
support the notion that the RWI is intrinsically two-dimensional. This implies that 
inconsistent upper disc boundary conditions, such as vanishing enthalpy perturbation, 
may inhibit the RWI in 3D. 

Il998l ) . The linear problem becomes a set of ordinary differen- 
tial equations (ODE), but these ODEs can be complicated. 
Furthermore, this approach does not p ermit boun dary con- 
ditions at other heights to be imposed (|Linll2012al . hereafter 



The RWI was originally described in two-dimensional 
(2D) discs. Recent works show that it can operate in three- 
dimensional (3 D) geometr y (jUm urhan 2010; Meh eut et al.l 
I2OIOI . l2012al lbl: iLinl l2012al ). In these calculations the per- 
turbation extends vertically through the disc, but the ef- 
fect of conditions at the disc surface was not investi- 
gated. This might be a relevant issue for the RWI to de- 
velop in non-magnetized regions of protoplanetary discs 
l|Lyra fc Mac Low! |2012| ). In the layered accretion sce- 
nario, this 'dead zone' is overlaid by an actively accreting 
layer (|Gammiell 19961 : lOkuzumi fc Hirose|[20Tll : lMartin et al.l 
|2012| ). A first step towards revealing vertical boundary ef- 
fects on the 3D RWI is to calculate it subject to different 
upper disc boundary conditions. 

Previous linear RWI simulations assume the zero- 
density surface as the upper disc boundary. When the den- 
sity stratification has the appropriate functional form, reg- 
ularity conditions enable the vertical dependence of pertur- 
bations t o be expressed in terrns of cl a ssic orthogonal poly- 
nomials (jPapaloizou fc Pringld 1 19851 : iTakeuchi fc Mivamal 



EH). 

Th e purpose of this work is to address the above caveat 
of IL12| . An alternative numerical method is applied to the 
linear problem so that effects of upper disc boundary condi- 
tions can be explored through simple numerical experiments. 
This paper is organized as follows. 32] defines the equilib- 
rium and perturbed disc. |3] gives a step-by-step description 
to convert the problem into a single matrix equation. Solu- 
tions are presented in 33] and discussed in ^S] 



2 MODEL AND GOVERNING EQUATIONS 

The system is a geometrically thin, non-self-gravitating, in- 
viscid fluid disc orbiting a central star of mass M«. The 
disc is governed by the Euler equations in (r, 2) cylindri- 
cal co-ordinates centered on the star. Units such that the 
gravitational constant G = M« = 1 are adopted. A poly- 
tropic equation of state is assumed, so that the pressure P 
is related to the density p by P = Kp^^^^^ , where n is the 

polytropic index and K \s & constant. 

Details of the equilibrium disc are given in 1L12| . The 
unperturbed disc is steady and axisymmetric with surface 
density profile in the form E oc r~°'B{r), where B{r) de- 
scribes a Gaussian bump centered at r = ro with ampli- 
tude A and width Ar. The density is p{r,z) = po('")[l — 
z'^ / [r]]'^ , where po is the midplane density and H{r) 
is the disc thickness. The unperturbed velocity field is 
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{vr,v^,Vz) = (0, rf2, 0) where ^{r) is the angular velocity, 
and Q, ~ Qk = \/ GM, j . Fiducial parameter values are: 
Q = 0.5, ro = 1, A = 1.4, Ar = 0.05ro and Ji'(ro) = 0.14ro. 

Eulerian perturbations to the above equilibrium in the 
form Re[5p(r, z) expi(crf + m(/>)] are considered, a = — (lj + 
17) is a complex eigenfrequency and m is a positive integer. 
7 and — tj are the mode growth rate and real frequency, re- 
spectively. The linearized fluid equations yield the following 
partial differential equation (PDE): 
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l|L12l ). where W = 5P/p is the enthalpy perturbation and 
D = — a^; where = r~^ d{r'^Q?) / dr is the square of the 
epicycle frequency and o" = a + mfi is the shifted frequency. 
The RWI is associ ated with an extremu m in potential vortic- 
ity (?7 = KV2nS. [Lovelace et al.lll999t) . In the above setup 
min(77) occurs near ro, so that a'(ro) ~ 0. The aim is to solve 
Eq. [T] as a PDE eigenvalue problem with various boundary 
conditions applied a,t z/H = Zg. 

For simplicity Zs is assumed constant, but extension 
to Zs = Zair) is straight forward provided that Zg(r) 7^ 
(in which case transformation to pola r-like coordinates is 
needed, see IPapaloizou fc Pringldfl^ '). 



2.1 Vertical boundary conditions 

The enthalpy perturbation is assumed to be symmetric 
about the midplane, since even m odes have been fou nd to 
dominate in nonlinear simulations (|Meheut et al.llioiol ) . The 
upper disc boundary is set to Zs < 1, thus avoiding the 
disc surface where the PDE becomes singular. The follow- 
ing upper boundary conditions (BC) were found to permit 
the RWI. 

'Open' boundary. The Lagrangian pressure perturbation 
vanishes at Zs'. 



(2) 



where ^ is the Lagrangian displacement. If Zs is unity then 
this is the usual condition for polytropes with vanishing 
pressure at its surface (which is satisfied automatically for 
regular solutions). 

'Solid' boundaries. At Zs the meridional velocity per- 
turbation satisfies 



(3) 



where u — i/{r) is a real function. Choosing 1/ — ZsdH/dr 
corresponds to zero velocity perpendicular to the upper disc 
boundary {Sv± — 0). This applies to an impermeable sur- 
face. The simpler choice = corresponds to zero vertical 
velocity {Sv^ — 0). This might apply, for example, when ver- 
tical velocities are strongly damped above Zs. In practice, 
the two conditions are similar because |!/| ^ 1 for a thin 
disc. 



3 NUMERICAL METHOD 

It is convenient to adopt the co-ordinate system {R, Z) = 
(r,z/H) so the computational domain is rectangular for con- 



stant Zs . Then the governing equation becomes 



A= (l-z2)i?^ 

B = -2Z{1 - Z^)R'^^, 

c.z',i-.',«'(f)'-i£igL_£! 



TTl 

Z'^)R^ + 2nZ^R^ — , 
H 



R^D 2 fH' 



1 - (2n+l)Z^ 



2mRn 



In 



2,. nPp^ R 



(4) 
(5) 

(6) 
(7) 
(8) 

(9) 

(10) 



and primes denote differentiation with respect to the argu- 
ment. 

To construct a numerical scheme, the radial co-ordinate 
is discretized into Nr grid points uniformly spaced by AR. 
Let Ri denote the radial co-ordinate of the i"^ grid point, 
define Wi{Z) as the solution along the vertical line R = 
Ri and approximate it with a linear combination of basis 
functions, 



Nz 



Wi{Z) = W{R^, Z) = ^afc,7/>fc(Z/Z. 



(11) 



for 1 ^ i ^ Nr, with = T2(fe_i), where T; is 
a Chebyshev polynomial of the first kind of order / 
(Abramowitz & Stogun 1965). The number of basis func- 
tions is Nz and /max = 2{Nz ~ 1) is the highest polynomial 
order. 

The pseudo-spectral coefficients Uki are obtained by de- 
manding the governing PDE, Eq.|3J to be satisfied at vertical 
grid points Zj along each line, where 



-Zs cos 



(j - 1 + «max/2)7r 



1 s: j s; iVz 



(12) 



This is a standard choice for collocation point s in pseudo - 
spectral methods with Chebyshev polynomials (|Bovdll200ll '). 
Note that only the upper half plane {z ^ 0) is covered, by 
the assumption of symmetry. 

The next step is to insert Eq.[TT]into Eq.|4]and evaluate 
the governing PDE at each {Ri,Zj). Radial derivatives are 
approximated with central differences while vertical deriva- 
tives are computed exactly. At each grid point the result 
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Table 1. Summary of linear simulations. The azimuthal 
wavenumber is m = 3 for all cases. 



(13) 



where Aij = A{Ri,Zj) and similarly for other PDE coef- 
ficients, and tl^kj = tl^kiZj /Zs). Eq. [13] can be written in 
matrix notation: 



0. 



(14) 



The operators V''' and are Nz x Nz matrices and the 
column vector a^'-' consists of the Nz pseudo-spectral coef- 
ficients associated with the i**^ line. V^*' represent coupling 
to adjacent lines. The matrix elements can be read off Eq. 

\m 

For simplicity, OrW — is set at radial boundariefl 
This is achieved implicitly by setting afc_i_i — ajc^i+i for 
i = 1, Nn. V^^-' and vi^^' are modified accordingly. At Zs, 
the governing equation is replaced by upper disc boundary 
conditions listed in ^2.1\ This corresponds to modifying the 
PDE coefficients and the matrix elements for j = Nz- 

There are Nh lines in total, so the entire system of 
equations is 



Ma = 0, 



(15) 



where M is a (NrNz) x (NrNz) block tridiagonal matrix 
and a is a column vector of length NrNz, representing the 
linear operators and the a^i, respectively. 

The discretized problem, Eg. 1151 is a set of linear homo- 
geneous equations. A solution method is described in IL12| . 
The eigenfrequency a is varied, using Newton- Raphson iter- 
ation, until the matrix M becomes singular. Only solutions 
where the reciprocal condition number of M is zero at ma- 
chine precision are accepted. LAPACK was used for the matrix 
operations, which also yield a up to an arbitrary normaliza- 
tion. 



Case 


n 




BC 




ui/mQo 






la 


1.5 


0.90 


AP = 





0.9941 


0.1074 


0.31 


lb 


1.5 


0.45 


AP = 





0.9997 


0.0889 


0.36 


2a 


1.5 


0.90 


5v±^ = 





0.9940 


0.1091 


0.26 


2b 


1.5 


0.45 


Svx^ = 





0.9939 


0.1174 


0.06 


3a 


1.5 


0.90 


Svz = 





0.9945 


0.1086 


0.26 


3b 


1.5 


0.45 


SVz = 





0.9973 


0.1143 


0.06 


4a 


0.7 


0.90 


AP = 





0.9959 


0.1976 


0.34 


4b 


0.7 


0.90 


SVz = 





0.9968 


0.2034 


0.22 


5a 


2.5 


0.90 


AP = 





0.9926 


0.0359 


0.39 


5b 


2.5 


0.90 


SVz = 





0.9927 


0.0361 


0.38 



{6,n), is a measure of three-dimensionality of the flow near 
ro. 6m is defined by 



(16) 



where Vr,z = Re[5vr,zW* {ro,0)]/\W{ro,0)\ is the real ve- 
locity perturbation at an azimuth close to a local maxi- 
mum in midplane enthalpy perturbation (hereafter the vor- 
tex core). {■) denotes spatial averaging over R £ [0.8, 1.2]ro 
and Z G [0, Z,]. 



4.1 Fiducial case 

The numerical method is first tested against IL12| with the 
open boundary (case la). Fig. [T] shows W at several heights 
and the perturbed real velocity field projected onto the (r, z) 
plane. \ W\ has a weak vertical dependence in the co-rotation 
region r £ [0.8, 1.2]ro. The enthalpy perturbation becomes 
increasingly three-dimensional away from ro. Near Zs, the 
disturbance in the wave-region can reach amplitudes compa- 
rable to that near ro. This means the RWI becomes radially 
global at large heights. 

The above r esult, considered as the intrinsic solution, 
is very similar to IL12| . including the growth rate and radi- 
ally converging flow towards ro with upwards motion at the 
vortex core. This validates the choice of solution method. 



4 LINEAR SIMULATIONS 

The computational domain is i? G [0.4, 1.6]ro and Zs = 0.9 
or Zs = 0.45. The default polytropic index is n = 1.5 but 
cases with n — 0.7, 2.5 are also presented. Linear modes 
with azimuthal wavenumber m — 3 are considered. 

The flducial case below employed a trial frequency 
a = -(0.994m + O.107i)r2o where fio = r2(ro). This was 
the solution found in ,L12 , Eigenfrequencies for other pa- 
rameter values may be found by varying those parame- 
ters slowly away from fiducial values. The numerical grid 
is (Nii,Nz ) = ( 384,8), which corresponds to imax = 14 (cf. 
^max = 6 in|L]J). 

Table [1] summarizes the calculations. The last column, 

^ The RWI typically has disturbance radially confined near tq 
and therefore insensitive to distant radial boundaries. 



4.2 Solid boundaries 

Comparison of case 2a ((5wx = 0) and 3a (Jwz = 0) with the 
reference case above show that a solid boundary only has a 
minor effect on the instability growth rate. The perturbation 
for case 3a is plotted in Fig. O 

In comparison with case la, the enthalpy perturbation 
in the co-rotation region is unaffected by solid boundaries. 
There is, however, a decrease in three-dimensionality as in- 
dicated by (^m). This correlates with a slight increase in 
growth rate, suggesting that suppressing vertical motions 
favor instability. 

A solid upper boundary causes the disturbance to be- 
come more two-dimensional in the wave-region. Compared 
to Fig. [TJ for r < 0.6ro the amplitude of VK(r, 0.9//) has de- 
creased while that for VF(r, 0) has increased. Consequently, 
the magnitude of the vertical flow is smaller in case 3a than 
in case la. 
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Figure 1. Reference case with AP = at Zs. The normalized en- 
thalpy perturbation, \W{r, z)/W{ro,0)\ (top), and the perturbed 
meridional velocity (bottom) are shown. The real velocity field is 
taken at the vortex core. Velocities are normalized so its magni- 
tude can be compared to other cases. 
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Figure 2. Same as Fig. [T] but with Svz = applied at Zs 




Figure 3. Midplane enthalpy perturbations subject to 5v^ 
at two values of Zs. 



4.3 Effect of domain height 

In case 2b and 3b the vertical domain is reduced from cases 
2a and 3a, respectively. This increases the growth rate by 
< 8% but {6m.) is reduced by a factor of 3 — 4. When 6m is av- 
eraged over Z £ [0, 0.45], cases 2a and 3a have {6m) ~ 0.19. 
So the reduction in three-dimensionality with decreased Zs 
is not simply due to averaging over a region closer to the 
midplane (towards which Svz — >■ 0). 

Fig. [3] compares |W^(r, 0)| between cases 2a and 2b. The 
co-rotation region is unchanged. (Note that Svz oc dzW/a 
and \a{ro)\/^lo ^ 1, so a small difference in dzW be- 
tween the two cases can still result in an appreciable dif- 
ference in vertical velocity.) The perturbation magnitude in 
r < 0.8, r > 1.2 has noticeably increased with reduced disc 
thickness. In this sense the RWI has been made more radi- 
ally global by restricting the fluid to a thinner slab. 

Case lb has AP = applied Zs = 0.45. This cor- 
responds to a situation where conditions above the upper 
boundary adjusts according to the RWI in Z < Zs, which 
is of course unrealistic. Nevertheless, it is interesting to 
note that the growth rate has decreased compared to case 
la, but as above this correlates with an increase in three- 
dimensionality. 



4.4 Dependence on polytropic index 



|L12| found that the vertical flow in the vortex core increases 
in magnitude with decreasing polytropic index when other 
parameters are fixed. This is reflected in Fig. U when the 
open boundary is applied. As n is lowered from 2.5 to 0.7, 
the maximum vertical velocity increases by a factor of ~ 6. 
For each n, a calculation with a solid boundary is also shown, 
for which the increase in vertical velocity is less significant 
(by a factor ~ 3). 

Fig.[4]indicate upper disc boundary conditions are more 
important for smaller polytropic indices. The dotted curve 
for n = 2.5 only deviates from the solid curve in Z > 0.6. 
For n — 0.7, forcing the flow to be horizontal at Zs causes 
the disturbance to be more two-dimensional throughout the 
column of the fluid. 

In previous studies with a vertically isothermal disc 
(equivalent to a p olytrope with n — >■ oo), the vortex core ap- 
pears hydrostatic (|Meheut et al.|[2012bl . lLl^ ). So for large n 
explicitly setting Sv^ = at makes little difference. Con- 
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Figure 4. Normalized vertical velocity near a local maximum in 
midplano enthalpy perturbation for a range of polytropic indices 
n and upper disc boundary conditions. Here v^. = rH^.. 



versely, as n is lowered the intrinsic solution acquires larger 
\vz\ at the vortex core, so that the effect of forcing Svz 
at the upper boundary becomes global in the vertical direc- 
tion. 

fmagine setting up the RWI with AP — a,t Za, 
but during its development force Sv^ — >■ at Z^. Sup- 
pose upper boundary effects are communicated through the 
fluid by sound waves. The sound-speed in a polytrope is 
Cs = ^JdP/dp = Q.kH{l - Z^)^''V\/2n, and the vertical 
sound-crossing time is Tc ~ H/cs oc y/n/Q.k- In unit time, 
the change in boundary condition should affect the lower n 
polytrope deeper into the disc than higher n, as observed in 
Fig- 131 (This is also seen when the bump amplitude A for 
each n was adjusted to achieve equal growth rates.) 



4.5 Other experiments 

Additional calculations have been performed with a vari- 
able surface function Zs = Zsir). Examples include a sur- 
face of constant aspect-ratio (with 5v±_ — 0) and a surface 
of constant pressure (with AP — 0). In this case a sec- 
ond co-ordinate transformation is required. Similar results 
to previous experiments were obtained, including the afore- 
mentioned trend of growth rate with respect to disc thick- 
ness. The RWI appears insensitive to the exact shape of the 
upper disc boundary as well. 



5 DISCUSSION 

Linear simulations of the RWI with different upper disc 
boundary conditions have been performed. Changing from 
a free surface to a solid boundary increase growth rates 
slightly. This is accompanied by a decrease in the three- 
dimensionality of the flow at the vortex core. Confining the 
fiuid to a smaller vertical domain has the s ame effect. Thi s 
suggests 3D effects are weakly stabilizing (|Li et al.ll2003l ). 
However, upper disc boundary conditions affect the vertical 
flow throughout the vertical extent of the vortex core. This 
effect is stronger for decreasing polytropic index n, which 
corresponds to decreasing compressibility. 

The above results support previous interpretations that 



the R WI is largely two-dimensional (jUmurhanl [2OI0I : iLh] 
l2012al ). The linear growth rate and general flow pattern is 
insensitive to details of the upper disc boundary conditions, 
provided it ts conststent with the 2D solution. A counter- 
example is vanishing enthalpy perturbation at the upper 
disc boundarjjf]. Indeed, in this case no RWI-like modes were 
found (i.e. modes with enthalpy disturbances radially con- 
fined to the bump). Inappropriate upper disc boundary con- 
ditions may suppress the RWI as see n in linear and nonlin ear 
simulations performed to date fe.g. iMeheut et al]l2012ah . 

Experimenting with different vertical boundary condi- 
tions, as done here, is only a crude way to mimic possible 
boundary effects in a protoplanetary disc. However, failure 
to find the RWI for some conditions imply that whether 
or not the RWI can develop in protoplanetary discs de- 
pends on the disc ver tical structure, which can be compli- 
cated (lTerauemll2008h . Calculating the RWI for protoplan- 
etary discs will ultimately require multi-layer fluid models 
(|Umurhanll2012l ). 

The solution method in this paper supersedes lL12l for its 
simplicity and ability to admit different upper disc boundary 
conditions. Also, it does not rely on the functional form of 
the background stratification. In principle, it can be applied 
to non-barotropic perturbations since inclusion of the energy 
equation only changes the coefficients of the governing PDE. 
Assuming the RWI is not modified significantly, then this 
solution method should be suitable for discs with ent ropy 
gradie nts. This will be the subject of a follow-up paper (jLinl 
l2012bl . submitted). 
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